Gam models: summary
Birthweight
2007-2020
m_BW_bam5_2 <- bam(BW ~ s(GA_weeks) + s(mat_age, k=7) + s(birth_Y_M_num) + s(month) + s(mean_ssep2) + parity_cat + sex + Urban3 + Language + mother_nationality_cat2, data=bevn_eco_in7)
summary(m_BW_bam5_2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## BW ~ s(GA_weeks) + s(mat_age, k = 7) + s(birth_Y_M_num) + s(month) +
## s(mean_ssep2) + parity_cat + sex + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3335.2704 0.8361 3989.044 < 2e-16 ***
## parity_cat(1,2] 127.4692 0.8344 152.762 < 2e-16 ***
## parity_cat(2,3] 173.6260 1.2772 135.939 < 2e-16 ***
## parity_cat(3,20] 199.6977 2.2149 90.160 < 2e-16 ***
## sex2 -147.0101 0.7458 -197.128 < 2e-16 ***
## Urban3 0.6608 0.7906 0.836 0.403
## LanguageFrench -48.0924 0.9182 -52.374 < 2e-16 ***
## LanguageItalian -63.1878 2.0493 -30.834 < 2e-16 ***
## mother_nationality_cat2Africa 11.0168 2.2735 4.846 1.26e-06 ***
## mother_nationality_cat2Asia -17.9253 2.0596 -8.703 < 2e-16 ***
## mother_nationality_cat2Europe 32.5679 0.8501 38.310 < 2e-16 ***
## mother_nationality_cat2Northern America 65.0093 4.8546 13.391 < 2e-16 ***
## mother_nationality_cat2Oceania 63.4678 11.1288 5.703 1.18e-08 ***
## mother_nationality_cat2Southern and Central America 47.0493 2.7991 16.809 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(GA_weeks) 8.826 8.988 73217.048 <2e-16 ***
## s(mat_age) 5.178 5.693 20.286 <2e-16 ***
## s(birth_Y_M_num) 8.285 8.859 39.009 <2e-16 ***
## s(month) 5.893 7.053 3.183 0.0022 **
## s(mean_ssep2) 7.634 8.534 43.693 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.401 Deviance explained = 40.1%
## fREML = 8.1205e+06 Scale est. = 1.5259e+05 n = 1099328 gam.check(m_BW_bam5_2) 
##
## Method: fREML Optimizer: perf newton
## full convergence after 11 iterations.
## Gradient range [-0.07345475,0.05455442]
## (score 8120451 & scale 152594.7).
## Hessian positive definite, eigenvalue range [1.608974,549654.6].
## Model rank = 56 / 56
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(GA_weeks) 9.00 8.83 0.98 0.095 .
## s(mat_age) 6.00 5.18 0.99 0.325
## s(birth_Y_M_num) 9.00 8.28 0.98 0.095 .
## s(month) 9.00 5.89 1.02 0.965
## s(mean_ssep2) 9.00 7.63 1.03 0.995
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 plot(m_BW_bam5_2, select=1, ylim=c(0, 4000),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam5_2)[1])
plot(m_BW_bam5_2, select=2, ylim=c(3200, 3400),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam5_2)[1])
plot(m_BW_bam5_2, select=3, ylim=c(3300, 3380),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam5_2)[1] , xlab = "1 : 2007-01, 25: 2009-01, 50: 2011-02, 75: 2013-03, 100: 2015-04, 125: 2017-05, 150: 2019-06")
plot(m_BW_bam5_2, select=4, ylim=c(3325, 3345),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam5_2)[1])
plot(m_BW_bam5_2, select=5, ylim=c(3200, 3400),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam5_2)[1])
Adjusting for the effect of the huge dataset (n=1,099,368)
Cohen’s d formula:
\(d=\frac{beta}{SD}\) and \(SD=\sqrt{n}*se\)
\(d=\frac{beta}{\sqrt{n}*se}\)
- beta: estimate
- se: standard error
- n: population/observation size
- SD: standard deviation
#getting the CI 95% and d (interpret only for the non-smooth parameters)
beta <- coef(m_BW_bam5_2)
Vb <- vcov(m_BW_bam5_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
n <- 1099328
d <- (abs(beta)/(sqrt(n)*se))
my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
round(head(my.ci,14), 2)## beta lci uci d
## (Intercept) 3335.27 3333.63 3336.91 3.80
## parity_cat(1,2] 127.47 125.83 129.10 0.15
## parity_cat(2,3] 173.63 171.12 176.13 0.13
## parity_cat(3,20] 199.70 195.36 204.04 0.09
## sex2 -147.01 -148.47 -145.55 0.19
## Urban3 0.66 -0.89 2.21 0.00
## LanguageFrench -48.09 -49.89 -46.29 0.05
## LanguageItalian -63.19 -67.20 -59.17 0.03
## mother_nationality_cat2Africa 11.02 6.56 15.47 0.00
## mother_nationality_cat2Asia -17.93 -21.96 -13.89 0.01
## mother_nationality_cat2Europe 32.57 30.90 34.23 0.04
## mother_nationality_cat2Northern America 65.01 55.49 74.52 0.01
## mother_nationality_cat2Oceania 63.47 41.66 85.28 0.01
## mother_nationality_cat2Southern and Central America 47.05 41.56 52.54 0.02Interpretation of m_BW_bam5_2 model :
Why choosing m_BW_bam5_2? It has the lowest AIC and also lower k value for maternal age, which seems more appropriate by looking at the graph.
- Non-smooth terms
Increasing parity increases birthweight compared to first births: from (126g, 129g) at second parity, (171g, 176g) at third parity, to (195g, 204g) for parities over 3.
Being born female decreases by (-148g, -146g) compared to males.
Urbanicity has no effect on birthweight.
Language: Compared to babies born in German(+Romansh)-speaking Switzerland, being born in the French-speaking part decreases BW by (-50g,-46g), while being born in the Italian-speaking part decreases BW by (-67g,-59g).
Maternal nationality: compared to babies born from Swiss mothers, babies born from Asian mothers are lighter (-22g, -14g). Babies born from other-nationalities mothers are heavier: (7g, 16g) for Africa, (31g, 34g) for other European countries, (55g, 75g) for Northern America, (42g, 86g) for Oceania, and (42g-53g) for Southern/Central America.
- Smooth terms
GA: Obviously, birthweight increases with increasing gestational age, with a plateau at around 800g for babies between 20 and 25 gestational weeks, and another plateau at around 4000g for babies above 45 weeks.
Maternal age: birthweight increases non-linearly with maternal age: for younger mothers (<25yo) it increases almost linearly from 3280g to 3350g, then declines very slightly until around 3340g when the mother is 25-40yo. Older women have heavier babies, increasing from 3340g to around 3360g for women 40-50yo (keep in mind we have excluded a few mothers above 50 - improbable).
Time (month+year combined): globally, birthweight is decreasing between 2007 and 2019 from 3350 to 3330 which is a minor change, and then increases from 2019 to reach 3350g. There are some other minor variations throughtout the whole time but this doesnt seem very relevant: only 20g variation.
Seasonal/month: minor sinusoidal seasonal variation between 3330-3340g, with lower birthweight in March and August. These difference dont seem much relevant either (10g variation will probably not impact the baby’s health).
Mean ssep: Birthweight increases by increasing ssep, linearly between ssep=25-45 approx (3250-3340g), and then there is some minor birthweight increase for higher ssep, until reaching a BW of around 3360g. Overall, a +100g increase is probably worht mentioning.
2008-2010
Because birth_Y_M_num and month variables correlated too much, moth variable was not kept in the model. Mean ssep was set as linear variables (when included as a smooth, the graph clearly showed a linear association).
m_BW_bam_2009_2 <- bam(BW ~ s(GA_weeks) + s(mat_age, k=7) + s(birth_Y_M_num) + mean_ssep2 + parity_cat + sex + Urban3 + Language + mother_nationality_cat2, data=bevn_eco_in7_2008_10)
summary(m_BW_bam_2009_2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## BW ~ s(GA_weeks) + s(mat_age, k = 7) + s(birth_Y_M_num) + mean_ssep2 +
## parity_cat + sex + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3278.4522 6.4015 512.135 < 2e-16 ***
## mean_ssep2 0.8785 0.1043 8.419 < 2e-16 ***
## parity_cat(1,2] 122.1534 1.8847 64.813 < 2e-16 ***
## parity_cat(2,3] 172.2947 2.8838 59.746 < 2e-16 ***
## parity_cat(3,20] 196.5618 4.9231 39.926 < 2e-16 ***
## sex2 -146.2210 1.6800 -87.037 < 2e-16 ***
## Urban3 0.8197 1.7447 0.470 0.63846
## LanguageFrench -45.9768 2.0608 -22.310 < 2e-16 ***
## LanguageItalian -65.0958 4.3465 -14.977 < 2e-16 ***
## mother_nationality_cat2Africa 41.2721 6.0377 6.836 8.18e-12 ***
## mother_nationality_cat2Asia -15.2504 4.9412 -3.086 0.00203 **
## mother_nationality_cat2Europe 36.4670 1.9396 18.801 < 2e-16 ***
## mother_nationality_cat2Northern America 82.8570 10.6237 7.799 6.25e-15 ***
## mother_nationality_cat2Oceania 54.3020 24.9696 2.175 0.02965 *
## mother_nationality_cat2Southern and Central America 60.6874 6.0057 10.105 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(GA_weeks) 8.246 8.808 14713.688 < 2e-16 ***
## s(mat_age) 3.134 3.771 7.669 1.06e-05 ***
## s(birth_Y_M_num) 2.838 3.533 14.969 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.397 Deviance explained = 39.7%
## fREML = 1.6251e+06 Scale est. = 1.548e+05 n = 219790 plot(m_BW_bam_2009_2, select=1, ylim=c(0, 4000),shade = TRUE, shade.col = "lightblue",shift = coef(m_BW_bam_2009_2)[1]) #GA
plot(m_BW_bam_2009_2, select=2, ylim=c(3200, 3350),shade = TRUE, shade.col = "lightblue",shift = coef(m_BW_bam_2009_2)[1]) #mat age
plot(m_BW_bam_2009_2, select=3, ylim=c(3250, 3350),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam_2009_2)[1], xlab = "13 : 2008-01, 20: 2008-08, 30: 2009-06, 40: 2010-04, 50: 2011-02") # time
gam.check(m_BW_bam_2009_2)
##
## Method: fREML Optimizer: perf newton
## full convergence after 9 iterations.
## Gradient range [-0.002922242,0.001533332]
## (score 1625083 & scale 154802.2).
## Hessian positive definite, eigenvalue range [0.2911485,109886].
## Model rank = 39 / 39
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(GA_weeks) 9.00 8.25 0.99 0.38
## s(mat_age) 6.00 3.13 0.98 0.08 .
## s(birth_Y_M_num) 9.00 2.84 0.99 0.36
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 concurvity(m_BW_bam_2009_2, full=TRUE) #values are now fine.## para s(GA_weeks) s(mat_age) s(birth_Y_M_num)
## worst 0.9828236 0.015759637 0.13707246 0.005557451
## observed 0.9828236 0.006364574 0.02758171 0.005155801
## estimate 0.9828236 0.009119694 0.08566593 0.003170851 # Summary of estimates, CI95% and Cohen's d
beta <- coef(m_BW_bam_2009_2)
Vb <- vcov(m_BW_bam_2009_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
n <- 219790
d <- (abs(beta)/(sqrt(n)*se))
my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
round(head(my.ci,15), 2) ## beta lci uci d
## (Intercept) 3278.45 3265.90 3291.00 1.09
## mean_ssep2 0.88 0.67 1.08 0.02
## parity_cat(1,2] 122.15 118.46 125.85 0.14
## parity_cat(2,3] 172.29 166.64 177.95 0.13
## parity_cat(3,20] 196.56 186.91 206.21 0.09
## sex2 -146.22 -149.51 -142.93 0.19
## Urban3 0.82 -2.60 4.24 0.00
## LanguageFrench -45.98 -50.02 -41.94 0.05
## LanguageItalian -65.10 -73.61 -56.58 0.03
## mother_nationality_cat2Africa 41.27 29.44 53.11 0.01
## mother_nationality_cat2Asia -15.25 -24.94 -5.57 0.01
## mother_nationality_cat2Europe 36.47 32.66 40.27 0.04
## mother_nationality_cat2Northern America 82.86 62.03 103.68 0.02
## mother_nationality_cat2Oceania 54.30 5.36 103.24 0.00
## mother_nationality_cat2Southern and Central America 60.69 48.92 72.46 0.022017-2020
We also did not use the month variable because of correlation between month and time (birth_Y_M_num) variables.
m_BW_bam_covid_2 <- bam(BW ~ s(GA_weeks) + s(mat_age, k=7) + s(birth_Y_M_num) + s(mean_ssep2) + parity_cat + sex + Urban3 + Language + mother_nationality_cat2, data=bevn_eco_in7_2017_20)
summary(m_BW_bam_covid_2)##
## Family: gaussian
## Link function: identity
##
## Formula:
## BW ~ s(GA_weeks) + s(mat_age, k = 7) + s(birth_Y_M_num) + s(mean_ssep2) +
## parity_cat + sex + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3341.78975 1.51250 2209.446 < 2e-16 ***
## parity_cat(1,2] 132.59338 1.50324 88.205 < 2e-16 ***
## parity_cat(2,3] 180.82725 2.29301 78.860 < 2e-16 ***
## parity_cat(3,20] 202.25754 4.00225 50.536 < 2e-16 ***
## sex2 -144.71388 1.34661 -107.465 < 2e-16 ***
## Urban3 -0.09571 1.42256 -0.067 0.94636
## LanguageFrench -51.18569 1.63285 -31.348 < 2e-16 ***
## LanguageItalian -62.48479 3.94863 -15.824 < 2e-16 ***
## mother_nationality_cat2Africa -6.51636 3.76791 -1.729 0.08373 .
## mother_nationality_cat2Asia -16.40941 3.54593 -4.628 3.70e-06 ***
## mother_nationality_cat2Europe 27.23301 1.52212 17.892 < 2e-16 ***
## mother_nationality_cat2Northern America 59.93225 9.02825 6.638 3.18e-11 ***
## mother_nationality_cat2Oceania 60.90901 22.28696 2.733 0.00628 **
## mother_nationality_cat2Southern and Central America 38.21351 5.36683 7.120 1.08e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(GA_weeks) 8.504 8.906 22216.026 <2e-16 ***
## s(mat_age) 4.935 5.500 23.982 <2e-16 ***
## s(birth_Y_M_num) 5.636 6.787 9.851 <2e-16 ***
## s(mean_ssep2) 5.067 6.235 15.020 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.399 Deviance explained = 39.9%
## fREML = 2.4557e+06 Scale est. = 1.5057e+05 n = 332752 plot(m_BW_bam_covid_2, select=1, ylim=c(0, 4000),shade = TRUE, shade.col = "lightblue",shift = coef(m_BW_bam_covid_2)[1])
plot(m_BW_bam_covid_2, select=2, ylim=c(3200, 3400),shade = TRUE, shade.col = "lightblue",shift = coef(m_BW_bam_covid_2)[1])
plot(m_BW_bam_covid_2, select=3, ylim=c(3310, 3360),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam_covid_2)[1],xlab = "120 : 2016-12, 130: 2017-10, 140: 2018-08, 150: 2019-06, 160: 2020-04")
plot(m_BW_bam_covid_2, select=4, ylim=c(3200, 3400),shade = TRUE, shade.col = "lightblue", shift = coef(m_BW_bam_covid_2)[1])
gam.check(m_BW_bam_covid_2)
##
## Method: fREML Optimizer: perf newton
## full convergence after 12 iterations.
## Gradient range [-0.02069349,0.01034373]
## (score 2455716 & scale 150568.5).
## Hessian positive definite, eigenvalue range [0.9600157,166367].
## Model rank = 47 / 47
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(GA_weeks) 9.00 8.50 1.00 0.39
## s(mat_age) 6.00 4.93 1.03 0.99
## s(birth_Y_M_num) 9.00 5.64 1.00 0.39
## s(mean_ssep2) 9.00 5.07 1.00 0.54 concurvity(m_BW_bam_covid_2, full=TRUE) #values are now fine.## para s(GA_weeks) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst 0.8024701 0.018706682 0.10865948 0.0012092151 0.17776696
## observed 0.8024701 0.007194495 0.06217053 0.0006631716 0.08814200
## estimate 0.8024701 0.010875530 0.06838436 0.0009362815 0.08021785 # Summary of estimates, CI95% and Cohen's d
beta <- coef(m_BW_bam_covid_2)
Vb <- vcov(m_BW_bam_covid_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
n <- 332752
my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))## Warning in cbind(beta, lci = beta - 1.96 * se, uci = beta + 1.96 * se, d): number of rows of result
## is not a multiple of vector length (arg 4) round(head(my.ci,14), 2)## beta lci uci d
## (Intercept) 3341.79 3338.83 3344.75 1.09
## parity_cat(1,2] 132.59 129.65 135.54 0.02
## parity_cat(2,3] 180.83 176.33 185.32 0.14
## parity_cat(3,20] 202.26 194.41 210.10 0.13
## sex2 -144.71 -147.35 -142.07 0.09
## Urban3 -0.10 -2.89 2.70 0.19
## LanguageFrench -51.19 -54.39 -47.98 0.00
## LanguageItalian -62.48 -70.23 -54.74 0.05
## mother_nationality_cat2Africa -6.52 -13.90 0.87 0.03
## mother_nationality_cat2Asia -16.41 -23.36 -9.46 0.01
## mother_nationality_cat2Europe 27.23 24.25 30.22 0.01
## mother_nationality_cat2Northern America 59.93 42.24 77.63 0.04
## mother_nationality_cat2Oceania 60.91 17.23 104.59 0.02
## mother_nationality_cat2Southern and Central America 38.21 27.69 48.73 0.00Stillbirth
2007-2020
m_SB_bam2 <- bam(stillbirth~s(GA_weeks)+ s(mat_age, k=7) + s(month) + birth_Y_M_num + mean_ssep2 + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in6)
summary(m_SB_bam2)##
## Family: binomial
## Link function: logit
##
## Formula:
## stillbirth ~ s(GA_weeks) + s(mat_age, k = 7) + s(month) + birth_Y_M_num +
## mean_ssep2 + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.619e+00 1.456e-01 -45.452 < 2e-16 ***
## birth_Y_M_num 4.545e-05 3.972e-04 0.114 0.90891
## mean_ssep2 -6.094e-03 2.332e-03 -2.614 0.00896 **
## Urban3 3.814e-02 3.892e-02 0.980 0.32712
## LanguageFrench 8.110e-02 4.404e-02 1.842 0.06554 .
## LanguageItalian -2.043e-02 1.024e-01 -0.200 0.84185
## mother_nationality_cat2Africa 2.275e-01 9.350e-02 2.434 0.01495 *
## mother_nationality_cat2Asia -3.683e-02 1.016e-01 -0.362 0.71708
## mother_nationality_cat2Europe 3.175e-02 4.253e-02 0.747 0.45532
## mother_nationality_cat2Northern America -3.488e-01 2.851e-01 -1.223 0.22121
## mother_nationality_cat2Oceania -1.733e-01 7.207e-01 -0.241 0.80993
## mother_nationality_cat2Southern and Central America -9.889e-02 1.351e-01 -0.732 0.46426
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 7.589 8.143 16734.305 <2e-16 ***
## s(mat_age) 3.147 3.820 13.061 0.0105 *
## s(month) 2.627 3.267 7.216 0.0701 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.133 Deviance explained = 32.1%
## fREML = 1.5655e+06 Scale est. = 1 n = 1102617 plot(m_SB_bam2, trans = plogis, select=1, shift = coef(m_SB_bam2)[1], shade = TRUE, shade.col = "lightblue") #GA
plot(m_SB_bam2, trans = plogis, select=2, ylim=c(0, 0.003), shift = coef(m_SB_bam2)[1], shade = TRUE, shade.col = "lightblue") #mat age
plot(m_SB_bam2, trans = plogis, select=3, ylim=c(0.0005, 0.0015), shift = coef(m_SB_bam2)[1], shade = TRUE, shade.col = "lightblue") #month
k.check(m_SB_bam2) ## k' edf k-index p-value
## s(GA_weeks) 9 7.588646 0.9509095 0.6175
## s(mat_age) 6 3.146678 0.9511237 0.6800
## s(month) 9 2.626666 0.9507136 0.5325 concurvity(m_SB_bam2)## para s(GA_weeks) s(mat_age) s(month)
## worst 0.9829412 0.010976422 0.05328595 0.0045643900
## observed 0.9829412 0.005611189 0.01081304 0.0005238286
## estimate 0.9829412 0.007092571 0.03420192 0.0034099344Calculation of Cohen’s d for OR:
\(d= \log(OR) * \frac{\sqrt{3}}{\pi}\)
\(d = \frac{log(OR)}{1.81}\)
#table to show OR, confidence interval and Cohen's d
OR <- exp(coef(m_SB_bam2))
beta <- coef(m_SB_bam2)
Vb <- vcov(m_SB_bam2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 12), 2)## OR lci uci d
## (Intercept) 0.00 0.00 0.00 3.66
## birth_Y_M_num 1.00 1.00 1.00 0.00
## mean_ssep2 0.99 0.99 1.00 0.00
## Urban3 1.04 0.96 1.12 0.02
## LanguageFrench 1.08 0.99 1.18 0.04
## LanguageItalian 0.98 0.80 1.20 0.01
## mother_nationality_cat2Africa 1.26 1.05 1.51 0.13
## mother_nationality_cat2Asia 0.96 0.79 1.18 0.02
## mother_nationality_cat2Europe 1.03 0.95 1.12 0.02
## mother_nationality_cat2Northern America 0.71 0.40 1.23 0.19
## mother_nationality_cat2Oceania 0.84 0.20 3.45 0.10
## mother_nationality_cat2Southern and Central America 0.91 0.70 1.18 0.05Time variables (birth_Y_M_num) and mean ssep variable were put as linear terms, because the model including them as smoothed terms showed that the relationship between these and stillbirth was linear.
#Adding BW –> discussing this.
m_SB_bam3 = update(m_SB_bam2, . ~ . + s(BW))
plot(m_SB_bam3, trans = plogis, select=4, ylim=c(0, 1), shift = coef(m_SB_bam3)[1], shade = TRUE, shade.col = "lightblue") #BW with whole CI
plot(m_SB_bam3, trans = plogis, select=4, ylim=c(0, 0.15), shift = coef(m_SB_bam3)[1], shade = TRUE, shade.col = "lightblue") #BW
concurvity(m_SB_bam3)## para s(GA_weeks) s(mat_age) s(month) s(BW)
## worst 0.9829455 0.8787415 0.05483844 0.0045794820 0.8787530
## observed 0.9829455 0.4617487 0.01746103 0.0006330721 0.5160650
## estimate 0.9829455 0.4140313 0.03533273 0.0034246658 0.3330819 k.check(m_SB_bam3)## k' edf k-index p-value
## s(GA_weeks) 9 7.293198 0.9949575 1.0000
## s(mat_age) 6 2.760090 0.9692146 0.5475
## s(month) 9 2.699310 0.9677931 0.2725
## s(BW) 9 6.933464 0.9951825 1.0000When we add birthweight, there is too much correlation betwwen gestational age and birthweight (concurvity shows 0.88 which is to close to 1): it is not a good idea to include BW in the model. But we can have a look at the smoothed curve for stillbirth depending on birthweight: it has a clear U shaped, with high stillbirth risk (2-11%) for very low birthweight babies (<2kg), and moderately high stillbirth risk (5-8%) for very high birthweight babies (>8kg).
2008-2010
m_SB_bam_2009_2<-bam(stillbirth~s(GA_weeks)+ s(mat_age, k=7) + birth_Y_M_num + mean_ssep2 + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in6_2008_10)
summary(m_SB_bam_2009_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## stillbirth ~ s(GA_weeks) + s(mat_age, k = 7) + birth_Y_M_num +
## mean_ssep2 + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.408576 0.339228 -18.892 <2e-16 ***
## birth_Y_M_num -0.003879 0.004052 -0.957 0.3385
## mean_ssep2 -0.007007 0.005230 -1.340 0.1803
## Urban3 0.174863 0.086722 2.016 0.0438 *
## LanguageFrench 0.033909 0.100425 0.338 0.7356
## LanguageItalian 0.043843 0.208702 0.210 0.8336
## mother_nationality_cat2Africa -0.307892 0.262929 -1.171 0.2416
## mother_nationality_cat2Asia -0.623295 0.292794 -2.129 0.0333 *
## mother_nationality_cat2Europe 0.065891 0.093736 0.703 0.4821
## mother_nationality_cat2Northern America -0.832096 0.740930 -1.123 0.2614
## mother_nationality_cat2Oceania -8.197695 52.128818 -0.157 0.8750
## mother_nationality_cat2Southern and Central America -0.318662 0.314474 -1.013 0.3109
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 5.937 6.945 3434.41 <2e-16 ***
## s(mat_age) 2.034 2.551 3.07 0.301
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.149 Deviance explained = 32.5%
## fREML = 3.1243e+05 Scale est. = 1 n = 220464 plot(m_SB_bam_2009_2, trans = plogis, select=1, ylim=c(0, 0.9) ,shift = coef(m_SB_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #GA
plot(m_SB_bam_2009_2, trans = plogis, select=2, ylim=c(0, 0.005), shift = coef(m_SB_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #mat age
k.check(m_SB_bam_2009_2)## k' edf k-index p-value
## s(GA_weeks) 9 5.936556 0.8528645 0.0150
## s(mat_age) 6 2.033857 0.9027108 0.2975 concurvity(m_SB_bam_2009_2)## para s(GA_weeks) s(mat_age)
## worst 0.9840551 0.010398851 0.061160016
## observed 0.9840551 0.005436465 0.005779778
## estimate 0.9840551 0.006773190 0.035691703 OR <- exp(coef(m_SB_bam_2009_2))
beta <- coef(m_SB_bam_2009_2)
Vb <- vcov(m_SB_bam_2009_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 12), 2)## OR lci uci d
## (Intercept) 0.00 0.00 0.000000e+00 3.54
## birth_Y_M_num 1.00 0.99 1.000000e+00 0.00
## mean_ssep2 0.99 0.98 1.000000e+00 0.00
## Urban3 1.19 1.00 1.410000e+00 0.10
## LanguageFrench 1.03 0.85 1.260000e+00 0.02
## LanguageItalian 1.04 0.69 1.570000e+00 0.02
## mother_nationality_cat2Africa 0.73 0.44 1.230000e+00 0.17
## mother_nationality_cat2Asia 0.54 0.30 9.500000e-01 0.34
## mother_nationality_cat2Europe 1.07 0.89 1.280000e+00 0.04
## mother_nationality_cat2Northern America 0.44 0.10 1.860000e+00 0.46
## mother_nationality_cat2Oceania 0.00 0.00 6.497288e+40 4.53
## mother_nationality_cat2Southern and Central America 0.73 0.39 1.350000e+00 0.18As for the previous model (including all timing 2007-2020), we kept time(birth_Y_M_num) and mean_ssep2 as linear terms.
2017-2020
m_SB_bam_covid_2 <-bam(stillbirth~s(GA_weeks)+ mat_age + birth_Y_M_num + mean_ssep2 + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in6_2017_20)
summary(m_SB_bam_covid_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## stillbirth ~ s(GA_weeks) + mat_age + birth_Y_M_num + mean_ssep2 +
## Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.794848 0.471639 -14.407 < 2e-16 ***
## mat_age 0.007085 0.006871 1.031 0.302447
## birth_Y_M_num -0.002495 0.002496 -1.000 0.317501
## mean_ssep2 -0.001567 0.004273 -0.367 0.713770
## Urban3 -0.121121 0.071465 -1.695 0.090107 .
## LanguageFrench 0.168834 0.078861 2.141 0.032282 *
## LanguageItalian -0.244986 0.230082 -1.065 0.286978
## mother_nationality_cat2Africa 0.552349 0.149243 3.701 0.000215 ***
## mother_nationality_cat2Asia 0.192085 0.167090 1.150 0.250312
## mother_nationality_cat2Europe 0.022826 0.078368 0.291 0.770843
## mother_nationality_cat2Northern America 0.274686 0.439897 0.624 0.532344
## mother_nationality_cat2Oceania -8.133629 46.445173 -0.175 0.860983
## mother_nationality_cat2Southern and Central America 0.044082 0.242353 0.182 0.855667
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 7.421 8.111 4940 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.127 Deviance explained = 33.1%
## fREML = 4.7434e+05 Scale est. = 1 n = 333721 plot(m_SB_bam_covid_2, trans = plogis, select=1, ylim=c(0, 1) ,shift = coef(m_SB_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #GA
OR <- exp(coef(m_SB_bam_covid_2))
beta <- coef(m_SB_bam_covid_2)
Vb <- vcov(m_SB_bam_covid_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 13), 2)## OR lci uci d
## (Intercept) 0.00 0.00 0.000000e+00 3.75
## mat_age 1.01 0.99 1.020000e+00 0.00
## birth_Y_M_num 1.00 0.99 1.000000e+00 0.00
## mean_ssep2 1.00 0.99 1.010000e+00 0.00
## Urban3 0.89 0.77 1.020000e+00 0.07
## LanguageFrench 1.18 1.01 1.380000e+00 0.09
## LanguageItalian 0.78 0.50 1.230000e+00 0.14
## mother_nationality_cat2Africa 1.74 1.30 2.330000e+00 0.31
## mother_nationality_cat2Asia 1.21 0.87 1.680000e+00 0.11
## mother_nationality_cat2Europe 1.02 0.88 1.190000e+00 0.01
## mother_nationality_cat2Northern America 1.32 0.56 3.120000e+00 0.15
## mother_nationality_cat2Oceania 0.00 0.00 1.005865e+36 4.49
## mother_nationality_cat2Southern and Central America 1.05 0.65 1.680000e+00 0.02 k.check(m_SB_bam_covid_2)## k' edf k-index p-value
## s(GA_weeks) 9 7.420665 0.8755089 0.02 concurvity(m_SB_bam_2009_2)## para s(GA_weeks) s(mat_age)
## worst 0.9840551 0.010398851 0.061160016
## observed 0.9840551 0.005436465 0.005779778
## estimate 0.9840551 0.006773190 0.035691703Why we chose this model: Maternal age had edf values very close to one: we put it as a linear term instead of smoothed variable. As for the previous model (including all timing 2007-2020), we kept time (birth_Y_M_num) and mean_ssep2 variables as linear terms.
Gestational age in days (continuous)
2007-2020
m_GA_bam_days <- bam(GA_days ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(month) + s(mean_ssep2) + parity_cat + sex + Urban3 + mother_nationality_cat2 + Language, data=bevn_eco_in7)
summary(m_GA_bam_days)##
## Family: gaussian
## Link function: identity
##
## Formula:
## GA_days ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(month) +
## s(mean_ssep2) + parity_cat + sex + Urban3 + mother_nationality_cat2 +
## Language
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 276.04111 0.01752 15751.707 < 2e-16 ***
## parity_cat(1,2] -1.99202 0.01754 -113.596 < 2e-16 ***
## parity_cat(2,3] -2.49376 0.02687 -92.814 < 2e-16 ***
## parity_cat(3,20] -2.56052 0.04653 -55.033 < 2e-16 ***
## sex2 1.81258 0.01581 114.644 < 2e-16 ***
## Urban3 0.09180 0.01657 5.539 3.04e-08 ***
## mother_nationality_cat2Africa 0.60971 0.04764 12.798 < 2e-16 ***
## mother_nationality_cat2Asia -0.79664 0.04316 -18.460 < 2e-16 ***
## mother_nationality_cat2Europe -0.35922 0.01783 -20.147 < 2e-16 ***
## mother_nationality_cat2Northern America -0.28620 0.10175 -2.813 0.00491 **
## mother_nationality_cat2Oceania -0.32052 0.23325 -1.374 0.16939
## mother_nationality_cat2Southern and Central America -1.59307 0.05865 -27.161 < 2e-16 ***
## LanguageFrench 0.10488 0.01928 5.441 5.29e-08 ***
## LanguageItalian -0.20128 0.04297 -4.684 2.81e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(BW) 8.977 9.000 127518.85 <2e-16 ***
## s(mat_age) 8.044 8.641 347.60 <2e-16 ***
## s(birth_Y_M_num) 8.689 8.972 115.30 <2e-16 ***
## s(month) 6.958 8.043 12.76 <2e-16 ***
## s(mean_ssep2) 7.867 8.668 59.48 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.515 Deviance explained = 51.5%
## fREML = 3.8714e+06 Scale est. = 67.031 n = 1099328 #plot(m_GA_bam_days, all.terms = TRUE,pages=1, shade = TRUE, shade.col = "lightblue")
plot(m_GA_bam_days, select=1, shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_days)[1]) #BW
plot(m_GA_bam_days, select=2, ylim=c(265, 280),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_days)[1]) #mat age
plot(m_GA_bam_days, select=3, ylim=c(274, 278),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_days)[1], xlab = "1 : 2007-01, 25: 2009-01, 50: 2011-02, 75: 2013-03, 100: 2015-04, 125: 2017-05, 150: 2019-06") #birth_Y_M_num
plot(m_GA_bam_days, select=4, ylim=c(275, 277),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_days)[1]) #month
plot(m_GA_bam_days, select=5, ylim=c(272, 280),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_days)[1]) #ssep
gam.check(m_GA_bam_days)
##
## Method: fREML Optimizer: perf newton
## full convergence after 13 iterations.
## Gradient range [-0.001200072,0.001194492]
## (score 3871436 & scale 67.03148).
## Hessian positive definite, eigenvalue range [1.310202,549654.5].
## Model rank = 59 / 59
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(BW) 9.00 8.98 1.00 0.500
## s(mat_age) 9.00 8.04 0.99 0.360
## s(birth_Y_M_num) 9.00 8.69 0.98 0.085 .
## s(month) 9.00 6.96 1.01 0.815
## s(mean_ssep2) 9.00 7.87 1.01 0.775
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 concurvity(m_GA_bam_days, full=TRUE)## para s(BW) s(mat_age) s(birth_Y_M_num) s(month) s(mean_ssep2)
## worst 0.8015268 0.04840718 0.12110433 0.03003088 0.028058758 0.17917903
## observed 0.8015268 0.01911969 0.07090416 0.01041702 0.001316444 0.04934839
## estimate 0.8015268 0.03133441 0.08072505 0.01100951 0.021033963 0.08954157 #estimate, 95%CI and d
beta <- coef(m_GA_bam_days)
Vb <- vcov(m_GA_bam_days, unconditional = TRUE)
se <- sqrt(diag(Vb))
n <- 1099328
d <- (abs(beta)/(sqrt(n)*se))
my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
round(head(my.ci,14), 2)## beta lci uci d
## (Intercept) 276.04 276.01 276.08 15.02
## parity_cat(1,2] -1.99 -2.03 -1.96 0.11
## parity_cat(2,3] -2.49 -2.55 -2.44 0.09
## parity_cat(3,20] -2.56 -2.65 -2.47 0.05
## sex2 1.81 1.78 1.84 0.11
## Urban3 0.09 0.06 0.12 0.01
## mother_nationality_cat2Africa 0.61 0.52 0.70 0.01
## mother_nationality_cat2Asia -0.80 -0.88 -0.71 0.02
## mother_nationality_cat2Europe -0.36 -0.39 -0.32 0.02
## mother_nationality_cat2Northern America -0.29 -0.49 -0.09 0.00
## mother_nationality_cat2Oceania -0.32 -0.78 0.14 0.00
## mother_nationality_cat2Southern and Central America -1.59 -1.71 -1.48 0.03
## LanguageFrench 0.10 0.07 0.14 0.01
## LanguageItalian -0.20 -0.29 -0.12 0.002008-2010
m_GA_bam_2009 <- bam(GA_days ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) + parity_cat + sex + Urban3 + mother_nationality_cat2 + Language, data=bevn_eco_in7_2008_10)
summary(m_GA_bam_2009)##
## Family: gaussian
## Link function: identity
##
## Formula:
## GA_days ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) +
## parity_cat + sex + Urban3 + mother_nationality_cat2 + Language
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 275.69818 0.03934 7008.938 < 2e-16 ***
## parity_cat(1,2] -1.85794 0.03988 -46.593 < 2e-16 ***
## parity_cat(2,3] -2.27191 0.06109 -37.189 < 2e-16 ***
## parity_cat(3,20] -2.34282 0.10408 -22.510 < 2e-16 ***
## sex2 1.81033 0.03583 50.528 < 2e-16 ***
## Urban3 0.16255 0.03745 4.340 1.42e-05 ***
## mother_nationality_cat2Africa -0.17976 0.12735 -1.412 0.15810
## mother_nationality_cat2Asia -0.87949 0.10418 -8.442 < 2e-16 ***
## mother_nationality_cat2Europe -0.29497 0.04100 -7.194 6.30e-13 ***
## mother_nationality_cat2Northern America -0.58820 0.22410 -2.625 0.00867 **
## mother_nationality_cat2Oceania -0.14096 0.52659 -0.268 0.78894
## mother_nationality_cat2Southern and Central America -1.51662 0.12663 -11.976 < 2e-16 ***
## LanguageFrench -0.08757 0.04398 -1.991 0.04649 *
## LanguageItalian -0.20814 0.09197 -2.263 0.02363 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(BW) 8.872 8.994 24566.37 <2e-16 ***
## s(mat_age) 7.177 7.939 66.93 <2e-16 ***
## s(birth_Y_M_num) 6.183 7.340 31.00 <2e-16 ***
## s(mean_ssep2) 4.779 5.919 20.22 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.505 Deviance explained = 50.5%
## fREML = 7.7699e+05 Scale est. = 68.838 n = 219790 #plot(m_GA_bam_2009, all.terms = TRUE,pages=1, shade = TRUE, shade.col = "lightblue")
plot(m_GA_bam_2009, select=1, shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_2009)[1]) #BW
plot(m_GA_bam_2009, select=2, ylim=c(265, 280),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_2009)[1]) #mat age
plot(m_GA_bam_2009, select=3, ylim=c(274, 278),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_2009)[1], xlab = "13 : 2008-01, 20: 2008-08, 30: 2009-06, 40: 2010-04, 50: 2011-02") #birth_Y_M_num
plot(m_GA_bam_2009, select=4, ylim=c(272, 280),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_2009)[1]) #ssep
gam.check(m_GA_bam_2009)
##
## Method: fREML Optimizer: perf newton
## full convergence after 12 iterations.
## Gradient range [-0.000335078,0.0003341467]
## (score 776990.4 & scale 68.83792).
## Hessian positive definite, eigenvalue range [0.9452316,109886].
## Model rank = 50 / 50
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(BW) 9.00 8.87 1.00 0.585
## s(mat_age) 9.00 7.18 0.96 0.005 **
## s(birth_Y_M_num) 9.00 6.18 1.01 0.710
## s(mean_ssep2) 9.00 4.78 1.00 0.475
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 concurvity(m_GA_bam_2009, full=TRUE)## para s(BW) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst 0.7978125 0.04851831 0.13569927 0.005392540 0.18029328
## observed 0.7978125 0.01941720 0.11222580 0.004359862 0.07127053
## estimate 0.7978125 0.03114362 0.08700719 0.002902380 0.08241982 #estimate, 95%CI and d
beta <- coef(m_GA_bam_2009)
Vb <- vcov(m_GA_bam_2009, unconditional = TRUE)
n <- 219790
se <- sqrt(diag(Vb))
d <- (abs(beta)/(sqrt(n)*se))
my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
round(head(my.ci,14), 2)## beta lci uci d
## (Intercept) 275.70 275.62 275.78 14.95
## parity_cat(1,2] -1.86 -1.94 -1.78 0.10
## parity_cat(2,3] -2.27 -2.39 -2.15 0.08
## parity_cat(3,20] -2.34 -2.55 -2.14 0.05
## sex2 1.81 1.74 1.88 0.11
## Urban3 0.16 0.09 0.24 0.01
## mother_nationality_cat2Africa -0.18 -0.43 0.07 0.00
## mother_nationality_cat2Asia -0.88 -1.08 -0.68 0.02
## mother_nationality_cat2Europe -0.29 -0.38 -0.21 0.02
## mother_nationality_cat2Northern America -0.59 -1.03 -0.15 0.01
## mother_nationality_cat2Oceania -0.14 -1.17 0.89 0.00
## mother_nationality_cat2Southern and Central America -1.52 -1.76 -1.27 0.03
## LanguageFrench -0.09 -0.17 0.00 0.00
## LanguageItalian -0.21 -0.39 -0.03 0.002017-2020
m_GA_bam_covid <- bam(GA_days ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) + parity_cat + sex + Urban3 + mother_nationality_cat2 + Language, data=bevn_eco_in7_2017_20)
summary(m_GA_bam_covid)##
## Family: gaussian
## Link function: identity
##
## Formula:
## GA_days ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) +
## parity_cat + sex + Urban3 + mother_nationality_cat2 + Language
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 276.47736 0.03150 8777.856 <2e-16 ***
## parity_cat(1,2] -2.04786 0.03143 -65.164 <2e-16 ***
## parity_cat(2,3] -2.69664 0.04796 -56.228 <2e-16 ***
## parity_cat(3,20] -2.85860 0.08356 -34.209 <2e-16 ***
## sex2 1.80880 0.02837 63.757 <2e-16 ***
## Urban3 0.05005 0.02973 1.683 0.0923 .
## mother_nationality_cat2Africa 0.91687 0.07848 11.682 <2e-16 ***
## mother_nationality_cat2Asia -0.94863 0.07386 -12.843 <2e-16 ***
## mother_nationality_cat2Europe -0.45937 0.03173 -14.477 <2e-16 ***
## mother_nationality_cat2Northern America -0.38787 0.18813 -2.062 0.0392 *
## mother_nationality_cat2Oceania -1.11572 0.46440 -2.402 0.0163 *
## mother_nationality_cat2Southern and Central America -1.84999 0.11179 -16.549 <2e-16 ***
## LanguageFrench 0.42080 0.03422 12.296 <2e-16 ***
## LanguageItalian -0.11568 0.08240 -1.404 0.1604
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(BW) 8.956 8.999 38956.550 < 2e-16 ***
## s(mat_age) 7.036 7.817 164.828 < 2e-16 ***
## s(birth_Y_M_num) 3.077 3.826 4.686 0.00137 **
## s(mean_ssep2) 7.577 8.501 15.089 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.517 Deviance explained = 51.7%
## fREML = 1.1677e+06 Scale est. = 65.374 n = 332752 #plot(m_GA_bam_covid, all.terms = TRUE,pages=1, shade = TRUE, shade.col = "lightblue")
plot(m_GA_bam_covid, select=1, shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_covid)[1]) #BW
plot(m_GA_bam_covid, select=2, ylim=c(265, 280),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_covid)[1]) #mat age
plot(m_GA_bam_covid, select=3, ylim=c(274, 278),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_covid)[1], xlab = "120 : 2016-12, 130: 2017-10, 140: 2018-08, 150: 2019-06, 160: 2020-04") #birth_Y_M_num
plot(m_GA_bam_covid, select=4, ylim=c(272, 280),shade = TRUE, shade.col = "lightblue", shift = coef(m_GA_bam_covid)[1]) #ssep
gam.check(m_GA_bam_covid)
##
## Method: fREML Optimizer: perf newton
## full convergence after 12 iterations.
## Gradient range [-0.0002086092,0.000205571]
## (score 1167708 & scale 65.37401).
## Hessian positive definite, eigenvalue range [0.6299583,166367].
## Model rank = 50 / 50
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(BW) 9.00 8.96 1.00 0.42
## s(mat_age) 9.00 7.04 1.01 0.82
## s(birth_Y_M_num) 9.00 3.08 1.00 0.50
## s(mean_ssep2) 9.00 7.58 1.01 0.70 concurvity(m_GA_bam_covid, full=TRUE)## para s(BW) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst 0.8021077 0.04825138 0.1070555 0.0012721330 0.17792321
## observed 0.8021077 0.01913284 0.0402550 0.0004357399 0.02758032
## estimate 0.8021077 0.02994665 0.0678886 0.0010076568 0.08028872 #estimate, 95%CI and d
beta <- coef(m_GA_bam_covid)
Vb <- vcov(m_GA_bam_covid, unconditional = TRUE)
n <- 332752
se <- sqrt(diag(Vb))
d <- (abs(beta)/(sqrt(n)*se))
my.ci <- data.frame(cbind(beta, lci = beta-1.96*se, uci = beta + 1.96*se, d))
round(head(my.ci,14), 2)## beta lci uci d
## (Intercept) 276.48 276.42 276.54 15.22
## parity_cat(1,2] -2.05 -2.11 -1.99 0.11
## parity_cat(2,3] -2.70 -2.79 -2.60 0.10
## parity_cat(3,20] -2.86 -3.02 -2.69 0.06
## sex2 1.81 1.75 1.86 0.11
## Urban3 0.05 -0.01 0.11 0.00
## mother_nationality_cat2Africa 0.92 0.76 1.07 0.02
## mother_nationality_cat2Asia -0.95 -1.09 -0.80 0.02
## mother_nationality_cat2Europe -0.46 -0.52 -0.40 0.03
## mother_nationality_cat2Northern America -0.39 -0.76 -0.02 0.00
## mother_nationality_cat2Oceania -1.12 -2.03 -0.21 0.00
## mother_nationality_cat2Southern and Central America -1.85 -2.07 -1.63 0.03
## LanguageFrench 0.42 0.35 0.49 0.02
## LanguageItalian -0.12 -0.28 0.05 0.00Pre-term birth: yes/no
2007-2020
m_PTB_bam<-bam(PTB~s(BW)+ s(mat_age) + s(birth_Y_M_num) + s(month) +s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7)
summary(m_PTB_bam)##
## Family: binomial
## Link function: logit
##
## Formula:
## PTB ~ s(BW) + s(mat_age) + s(birth_Y_M_num) + s(month) + s(mean_ssep2) +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.27613 0.01550 -275.892 < 2e-16 ***
## parity_cat(1,2] 0.09992 0.01228 8.135 4.11e-16 ***
## parity_cat(2,3] 0.22338 0.01922 11.620 < 2e-16 ***
## parity_cat(3,20] 0.39947 0.03191 12.520 < 2e-16 ***
## Urban3 -0.05225 0.01142 -4.574 4.78e-06 ***
## LanguageFrench -0.15425 0.01311 -11.770 < 2e-16 ***
## LanguageItalian -0.19535 0.02835 -6.890 5.57e-12 ***
## mother_nationality_cat2Africa -0.13392 0.03577 -3.743 0.000182 ***
## mother_nationality_cat2Asia -0.10938 0.02933 -3.729 0.000192 ***
## mother_nationality_cat2Europe 0.05210 0.01243 4.192 2.77e-05 ***
## mother_nationality_cat2Northern America 0.12821 0.07285 1.760 0.078403 .
## mother_nationality_cat2Oceania 0.25209 0.15753 1.600 0.109543
## mother_nationality_cat2Southern and Central America 0.21445 0.03881 5.525 3.29e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(BW) 8.345 8.863 97491.56 < 2e-16 ***
## s(mat_age) 6.424 7.331 146.40 < 2e-16 ***
## s(birth_Y_M_num) 3.072 3.822 80.01 < 2e-16 ***
## s(month) 5.348 6.475 15.94 0.01554 *
## s(mean_ssep2) 5.483 6.674 23.79 0.00097 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.373 Deviance explained = 42.5%
## fREML = 1.5599e+06 Scale est. = 1 n = 1099328 #adding intercept
#plot(m_PTB_bam, pages = 1, trans = plogis, shift = coef(m_PTB_bam)[1])
plot(m_PTB_bam, trans = plogis, select=1,,shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue") #BW
plot(m_PTB_bam, trans = plogis, select=2, ylim=c(0, 0.05), shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue") #matage
plot(m_PTB_bam, trans = plogis, select=3, ylim=c(0.01, 0.02), shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue", xlab = "1 : 2007-01, 25: 2009-01, 50: 2011-02, 75: 2013-03, 100: 2015-04, 125: 2017-05, 150: 2019-06") #birth_Y_M_num
plot(m_PTB_bam, trans = plogis, select=4, ylim=c(0.012, 0.016), shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue") #month
plot(m_PTB_bam, trans = plogis, select=5, ylim=c(0.001, 0.02),shift = coef(m_PTB_bam)[1], shade = TRUE, shade.col = "lightblue") #ssep
k.check(m_PTB_bam) #some kindex values are a bit far from 1## k' edf k-index p-value
## s(BW) 9 8.345454 0.9139217 0.0525
## s(mat_age) 9 6.424215 0.9353811 0.4775
## s(birth_Y_M_num) 9 3.072089 0.9373408 0.5425
## s(month) 9 5.347513 0.9337870 0.4225
## s(mean_ssep2) 9 5.482747 0.9179691 0.1225 concurvity(m_PTB_bam)## para s(BW) s(mat_age) s(birth_Y_M_num) s(month) s(mean_ssep2)
## worst 0.75695 0.02679970 0.12110615 0.03003623 0.028060770 0.17918840
## observed 0.75695 0.02174550 0.09399625 0.01560116 0.001701081 0.06553150
## estimate 0.75695 0.01814849 0.08072682 0.01101228 0.021035075 0.08951894 #to calculate the OR of PTB
OR <- exp(coef(m_PTB_bam))
beta <- coef(m_PTB_bam)
Vb <- vcov(m_PTB_bam, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 13), 2)## OR lci uci d
## (Intercept) 0.01 0.01 0.01 2.36
## parity_cat(1,2] 1.11 1.08 1.13 0.06
## parity_cat(2,3] 1.25 1.20 1.30 0.12
## parity_cat(3,20] 1.49 1.40 1.59 0.22
## Urban3 0.95 0.93 0.97 0.03
## LanguageFrench 0.86 0.84 0.88 0.09
## LanguageItalian 0.82 0.78 0.87 0.11
## mother_nationality_cat2Africa 0.87 0.82 0.94 0.07
## mother_nationality_cat2Asia 0.90 0.85 0.95 0.06
## mother_nationality_cat2Europe 1.05 1.03 1.08 0.03
## mother_nationality_cat2Northern America 1.14 0.99 1.31 0.07
## mother_nationality_cat2Oceania 1.29 0.94 1.75 0.14
## mother_nationality_cat2Southern and Central America 1.24 1.15 1.34 0.122008-2010
The variables ssep, month and birth_Y_M_num were set as non-smooth variables (model including them as smooth variables showed a linear relationship with PTB risk)
m_PTB_bam_2009_2 <-bam(PTB~s(BW)+ s(mat_age) + birth_Y_M_num + month + mean_ssep2 + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2008_10)
summary(m_PTB_bam_2009_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## PTB ~ s(BW) + s(mat_age) + birth_Y_M_num + month + mean_ssep2 +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.046320 0.098396 -41.123 < 2e-16 ***
## birth_Y_M_num -0.005923 0.001215 -4.876 1.08e-06 ***
## month 0.002883 0.003648 0.790 0.429447
## mean_ssep2 0.001002 0.001454 0.689 0.490710
## parity_cat(1,2] 0.041571 0.027001 1.540 0.123654
## parity_cat(2,3] 0.180608 0.042228 4.277 1.89e-05 ***
## parity_cat(3,20] 0.338843 0.068907 4.917 8.77e-07 ***
## Urban3 -0.062921 0.024527 -2.565 0.010307 *
## LanguageFrench -0.124757 0.028516 -4.375 1.21e-05 ***
## LanguageItalian -0.233616 0.058557 -3.990 6.62e-05 ***
## mother_nationality_cat2Africa -0.016162 0.088601 -0.182 0.855263
## mother_nationality_cat2Asia -0.075412 0.068437 -1.102 0.270497
## mother_nationality_cat2Europe 0.081565 0.027526 2.963 0.003045 **
## mother_nationality_cat2Northern America 0.330363 0.148260 2.228 0.025863 *
## mother_nationality_cat2Oceania 0.481704 0.313962 1.534 0.124962
## mother_nationality_cat2Southern and Central America 0.273003 0.081406 3.354 0.000798 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(BW) 6.840 7.709 19924.46 <2e-16 ***
## s(mat_age) 5.576 6.569 42.04 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.362 Deviance explained = 41.2%
## fREML = 3.1193e+05 Scale est. = 1 n = 219790 plot(m_PTB_bam_2009_2, trans = plogis, select=1, shade = TRUE, shade.col = "lightblue", shift = coef(m_PTB_bam_2009_2)[1]) #BW
plot(m_PTB_bam_2009_2, trans = plogis, select=2, ylim=c(0, 0.1),shade = TRUE, shade.col = "lightblue", shift = coef(m_PTB_bam_2009_2)[1]) #mat age
k.check(m_PTB_bam_2009_2)## k' edf k-index p-value
## s(BW) 9 6.839574 0.9350836 0.0325
## s(mat_age) 9 5.576486 0.9457077 0.1000 concurvity(m_PTB_bam_2009_2)## para s(BW) s(mat_age)
## worst 0.9850434 0.02716942 0.13455667
## observed 0.9850434 0.02211902 0.09906189
## estimate 0.9850434 0.01766117 0.08612003 #to calculate the OR of PTB
OR <- exp(coef(m_PTB_bam_2009_2))
beta <- coef(m_PTB_bam_2009_2)
Vb <- vcov(m_PTB_bam_2009_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 16), 2)## OR lci uci d
## (Intercept) 0.02 0.01 0.02 2.24
## birth_Y_M_num 0.99 0.99 1.00 0.00
## month 1.00 1.00 1.01 0.00
## mean_ssep2 1.00 1.00 1.00 0.00
## parity_cat(1,2] 1.04 0.99 1.10 0.02
## parity_cat(2,3] 1.20 1.10 1.30 0.10
## parity_cat(3,20] 1.40 1.23 1.61 0.19
## Urban3 0.94 0.89 0.99 0.03
## LanguageFrench 0.88 0.83 0.93 0.07
## LanguageItalian 0.79 0.71 0.89 0.13
## mother_nationality_cat2Africa 0.98 0.83 1.17 0.01
## mother_nationality_cat2Asia 0.93 0.81 1.06 0.04
## mother_nationality_cat2Europe 1.08 1.03 1.15 0.05
## mother_nationality_cat2Northern America 1.39 1.04 1.86 0.18
## mother_nationality_cat2Oceania 1.62 0.87 3.00 0.27
## mother_nationality_cat2Southern and Central America 1.31 1.12 1.54 0.152017-2020
#birth_Y_M, month and mean ssep as linear variables:
m_PTB_bam_covid_2 <-bam(PTB~s(BW)+ s(mat_age) + birth_Y_M_num + month + mean_ssep2 + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2017_20)
summary(m_PTB_bam_covid_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## PTB ~ s(BW) + s(mat_age) + birth_Y_M_num + month + mean_ssep2 +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.3365644 0.1328509 -32.642 < 2e-16 ***
## birth_Y_M_num -0.0009770 0.0007532 -1.297 0.194627
## month 0.0008431 0.0030337 0.278 0.781074
## mean_ssep2 0.0006166 0.0012510 0.493 0.622082
## parity_cat(1,2] 0.1638479 0.0228156 7.181 6.9e-13 ***
## parity_cat(2,3] 0.2983101 0.0355576 8.389 < 2e-16 ***
## parity_cat(3,20] 0.5173515 0.0576642 8.972 < 2e-16 ***
## Urban3 -0.0637814 0.0208789 -3.055 0.002252 **
## LanguageFrench -0.2005760 0.0240309 -8.347 < 2e-16 ***
## LanguageItalian -0.1866281 0.0561558 -3.323 0.000889 ***
## mother_nationality_cat2Africa -0.2119198 0.0622213 -3.406 0.000659 ***
## mother_nationality_cat2Asia -0.0952343 0.0519765 -1.832 0.066913 .
## mother_nationality_cat2Europe 0.0586032 0.0228280 2.567 0.010254 *
## mother_nationality_cat2Northern America 0.0482582 0.1448537 0.333 0.739020
## mother_nationality_cat2Oceania 0.1284141 0.3413723 0.376 0.706790
## mother_nationality_cat2Southern and Central America 0.2060507 0.0754333 2.732 0.006303 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(BW) 8.033 8.714 28326.69 < 2e-16 ***
## s(mat_age) 4.845 5.832 35.55 3.13e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.377 Deviance explained = 43.3%
## fREML = 4.7393e+05 Scale est. = 1 n = 332752 plot(m_PTB_bam_covid_2, trans = plogis, select=1, shade = TRUE, shade.col = "lightblue", shift = coef(m_PTB_bam_covid_2)[1]) #BW
plot(m_PTB_bam_covid_2, trans = plogis, select=2, ylim=c(0, 0.05),shade = TRUE, shade.col = "lightblue", shift = coef(m_PTB_bam_covid_2)[1]) #mat age
k.check(m_PTB_bam_covid_2)## k' edf k-index p-value
## s(BW) 9 8.033448 0.9925878 0.9975
## s(mat_age) 9 4.844914 0.9445570 0.2325 concurvity(m_PTB_bam_covid_2)## para s(BW) s(mat_age)
## worst 0.9940987 0.02711559 0.10614541
## observed 0.9940987 0.02201811 0.04842693
## estimate 0.9940987 0.01778896 0.06719466 #to calculate the OR of PTB
OR <- exp(coef(m_PTB_bam_covid_2))
beta <- coef(m_PTB_bam_covid_2)
Vb <- vcov(m_PTB_bam_covid_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 16), 2)## OR lci uci d
## (Intercept) 0.01 0.01 0.02 2.40
## birth_Y_M_num 1.00 1.00 1.00 0.00
## month 1.00 0.99 1.01 0.00
## mean_ssep2 1.00 1.00 1.00 0.00
## parity_cat(1,2] 1.18 1.13 1.23 0.09
## parity_cat(2,3] 1.35 1.26 1.44 0.16
## parity_cat(3,20] 1.68 1.50 1.88 0.29
## Urban3 0.94 0.90 0.98 0.04
## LanguageFrench 0.82 0.78 0.86 0.11
## LanguageItalian 0.83 0.74 0.93 0.10
## mother_nationality_cat2Africa 0.81 0.72 0.91 0.12
## mother_nationality_cat2Asia 0.91 0.82 1.01 0.05
## mother_nationality_cat2Europe 1.06 1.01 1.11 0.03
## mother_nationality_cat2Northern America 1.05 0.79 1.39 0.03
## mother_nationality_cat2Oceania 1.14 0.58 2.22 0.07
## mother_nationality_cat2Southern and Central America 1.23 1.06 1.42 0.11Low birth weight
2007-2020
Month variable was set as a linear term (when included as a smooth term, edf=1.3)
m_LBW_bam<-bam(LBW~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) + month +s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7)
summary(m_LBW_bam)##
## Family: binomial
## Link function: logit
##
## Formula:
## LBW ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + month + s(mean_ssep2) +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.1815402 0.0201764 -207.249 < 2e-16 ***
## month 0.0006484 0.0017701 0.366 0.7141
## parity_cat(1,2] -0.6579649 0.0139870 -47.041 < 2e-16 ***
## parity_cat(2,3] -0.7937344 0.0226928 -34.977 < 2e-16 ***
## parity_cat(3,20] -0.8078470 0.0386216 -20.917 < 2e-16 ***
## Urban3 0.0172269 0.0125976 1.367 0.1715
## LanguageFrench 0.2850660 0.0140620 20.272 < 2e-16 ***
## LanguageItalian 0.1779457 0.0302849 5.876 4.21e-09 ***
## mother_nationality_cat2Africa 0.0007782 0.0376698 0.021 0.9835
## mother_nationality_cat2Asia -0.0152786 0.0317249 -0.482 0.6301
## mother_nationality_cat2Europe -0.1505381 0.0139476 -10.793 < 2e-16 ***
## mother_nationality_cat2Northern America -0.4943755 0.0911825 -5.422 5.90e-08 ***
## mother_nationality_cat2Oceania -0.4153060 0.2025663 -2.050 0.0403 *
## mother_nationality_cat2Southern and Central America -0.3650949 0.0457429 -7.981 1.45e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 7.933 8.460 89283.36 <2e-16 ***
## s(mat_age) 5.676 6.654 115.71 <2e-16 ***
## s(birth_Y_M_num) 3.911 4.836 55.59 <2e-16 ***
## s(mean_ssep2) 2.274 2.919 130.07 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.409 Deviance explained = 45.3%
## fREML = 1.5977e+06 Scale est. = 1 n = 1099328 plot(m_LBW_bam, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_bam)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks
plot(m_LBW_bam, trans = plogis, select=2, ylim=c(.005, .03), shift = coef(m_LBW_bam)[1], shade = TRUE, shade.col = "lightblue") #mat_age
plot(m_LBW_bam, trans = plogis, select=3, ylim=c(.01, .02), shift = coef(m_LBW_bam)[1], shade = TRUE, shade.col = "lightblue", xlab = "1 : 2007-01, 25: 2009-01, 50: 2011-02, 75: 2013-03, 100: 2015-04, 125: 2017-05, 150: 2019-06") #birth_Y_M_num
plot(m_LBW_bam, trans = plogis, select=4, ylim=c(.005, .025), shift = coef(m_LBW_bam)[1], shade = TRUE, shade.col = "lightblue") #ssep
k.check(m_LBW_bam)## k' edf k-index p-value
## s(GA_weeks) 9 7.933245 0.8882932 0.0000
## s(mat_age) 9 5.675613 0.9391889 0.1550
## s(birth_Y_M_num) 9 3.910952 0.9772231 0.9800
## s(mean_ssep2) 9 2.273894 0.9575589 0.5925 concurvity(m_LBW_bam)## para s(GA_weeks) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst 0.8741049 0.016894030 0.12294860 0.03012542 0.17905289
## observed 0.8741049 0.007885243 0.04784953 0.01087336 0.12485367
## estimate 0.8741049 0.009156089 0.08218594 0.01143670 0.08949728 #to calculate the OR of LBW for non smooth terms
OR <- exp(coef(m_LBW_bam))
beta <- coef(m_LBW_bam)
Vb <- vcov(m_LBW_bam, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 14), 2)## OR lci uci d
## (Intercept) 0.02 0.01 0.02 2.31
## month 1.00 1.00 1.00 0.00
## parity_cat(1,2] 0.52 0.50 0.53 0.36
## parity_cat(2,3] 0.45 0.43 0.47 0.44
## parity_cat(3,20] 0.45 0.41 0.48 0.45
## Urban3 1.02 0.99 1.04 0.01
## LanguageFrench 1.33 1.29 1.37 0.16
## LanguageItalian 1.19 1.13 1.27 0.10
## mother_nationality_cat2Africa 1.00 0.93 1.08 0.00
## mother_nationality_cat2Asia 0.98 0.93 1.05 0.01
## mother_nationality_cat2Europe 0.86 0.84 0.88 0.08
## mother_nationality_cat2Northern America 0.61 0.51 0.73 0.27
## mother_nationality_cat2Oceania 0.66 0.44 0.98 0.23
## mother_nationality_cat2Southern and Central America 0.69 0.63 0.76 0.202008-2010
As the variables birth_Y_M_num and month are too correlated : I did the model without month variable. Also, birth_Y_M_num and mean_ssep2 were set as non smooth variables (edf ~1 fit hey were included as smooth variables)
m_LBW_bam_2009_2 <-bam(LBW~s(GA_weeks)+ s(mat_age) + birth_Y_M_num + mean_ssep2 + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2008_10)
summary(m_LBW_bam_2009_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## LBW ~ s(GA_weeks) + s(mat_age) + birth_Y_M_num + mean_ssep2 +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.775443 0.109629 -34.438 < 2e-16 ***
## birth_Y_M_num 0.007496 0.001295 5.791 7.00e-09 ***
## mean_ssep2 -0.009640 0.001647 -5.853 4.82e-09 ***
## parity_cat(1,2] -0.630716 0.031012 -20.338 < 2e-16 ***
## parity_cat(2,3] -0.733919 0.050453 -14.547 < 2e-16 ***
## parity_cat(3,20] -0.783168 0.085537 -9.156 < 2e-16 ***
## Urban3 0.052812 0.027693 1.907 0.056512 .
## LanguageFrench 0.246728 0.031327 7.876 3.38e-15 ***
## LanguageItalian 0.296893 0.061089 4.860 1.17e-06 ***
## mother_nationality_cat2Africa -0.074333 0.095410 -0.779 0.435928
## mother_nationality_cat2Asia -0.024941 0.075948 -0.328 0.742615
## mother_nationality_cat2Europe -0.099987 0.031033 -3.222 0.001273 **
## mother_nationality_cat2Northern America -0.725315 0.208002 -3.487 0.000488 ***
## mother_nationality_cat2Oceania -0.138900 0.388893 -0.357 0.720967
## mother_nationality_cat2Southern and Central America -0.423089 0.099577 -4.249 2.15e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 7.110 7.829 18402.4 < 2e-16 ***
## s(mat_age) 3.831 4.745 21.7 0.00047 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.401 Deviance explained = 44.1%
## fREML = 3.1231e+05 Scale est. = 1 n = 219790 plot(m_LBW_bam_2009_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks
plot(m_LBW_bam_2009_2, trans = plogis, select=2, ylim=c(.005, .045), shift = coef(m_LBW_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age
k.check(m_LBW_bam_2009_2)## k' edf k-index p-value
## s(GA_weeks) 9 7.109732 0.9524535 0.5325
## s(mat_age) 9 3.831235 0.9455602 0.3725 concurvity(m_LBW_bam_2009_2)## para s(GA_weeks) s(mat_age)
## worst 0.984819 0.015475405 0.13702511
## observed 0.984819 0.007432975 0.03986371
## estimate 0.984819 0.008446321 0.08791625 #to calculate the OR of LBW for non smooth terms
OR <- exp(coef(m_LBW_bam_2009_2))
beta <- coef(m_LBW_bam_2009_2)
Vb <- vcov(m_LBW_bam_2009_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 15), 2)## OR lci uci d
## (Intercept) 0.02 0.02 0.03 2.09
## birth_Y_M_num 1.01 1.00 1.01 0.00
## mean_ssep2 0.99 0.99 0.99 0.01
## parity_cat(1,2] 0.53 0.50 0.57 0.35
## parity_cat(2,3] 0.48 0.43 0.53 0.41
## parity_cat(3,20] 0.46 0.39 0.54 0.43
## Urban3 1.05 1.00 1.11 0.03
## LanguageFrench 1.28 1.20 1.36 0.14
## LanguageItalian 1.35 1.19 1.52 0.16
## mother_nationality_cat2Africa 0.93 0.77 1.12 0.04
## mother_nationality_cat2Asia 0.98 0.84 1.13 0.01
## mother_nationality_cat2Europe 0.90 0.85 0.96 0.06
## mother_nationality_cat2Northern America 0.48 0.32 0.73 0.40
## mother_nationality_cat2Oceania 0.87 0.41 1.87 0.08
## mother_nationality_cat2Southern and Central America 0.66 0.54 0.80 0.232017-2020
Month variable was also not included in the model (correlation between birth_Y_M and month variables), and we set k=3 for s(birth_Y_M_num) because its almost linear (edf=1.99). Mean ssep set as linear term because almost linear (edf=1.47).
m_LBW_bam_covid_2 <-bam(LBW~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num, k=3) + mean_ssep2 + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2017_20)
summary(m_LBW_bam_covid_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## LBW ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num, k = 3) + mean_ssep2 +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.789475 0.086217 -43.953 < 2e-16 ***
## mean_ssep2 -0.008699 0.001378 -6.312 2.75e-10 ***
## parity_cat(1,2] -0.692985 0.025756 -26.905 < 2e-16 ***
## parity_cat(2,3] -0.864438 0.041814 -20.674 < 2e-16 ***
## parity_cat(3,20] -0.804935 0.068636 -11.728 < 2e-16 ***
## Urban3 0.021859 0.023021 0.950 0.34236
## LanguageFrench 0.293669 0.025656 11.446 < 2e-16 ***
## LanguageItalian 0.163459 0.059747 2.736 0.00622 **
## mother_nationality_cat2Africa 0.072592 0.063694 1.140 0.25441
## mother_nationality_cat2Asia -0.132195 0.057206 -2.311 0.02084 *
## mother_nationality_cat2Europe -0.187461 0.025499 -7.352 1.96e-13 ***
## mother_nationality_cat2Northern America -0.533030 0.172905 -3.083 0.00205 **
## mother_nationality_cat2Oceania -1.016028 0.490264 -2.072 0.03823 *
## mother_nationality_cat2Southern and Central America -0.263475 0.085099 -3.096 0.00196 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 7.293 7.922 26160.737 <2e-16 ***
## s(mat_age) 4.829 5.820 55.520 <2e-16 ***
## s(birth_Y_M_num) 1.705 1.913 2.865 0.163
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.41 Deviance explained = 46.1%
## fREML = 4.7418e+05 Scale est. = 1 n = 332752 plot(m_LBW_bam_covid_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks
plot(m_LBW_bam_covid_2, trans = plogis, select=2, ylim=c(.005, .06), shift = coef(m_LBW_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age
plot(m_LBW_bam_covid_2, trans = plogis, select=3, ylim=c(.015, .025), shift = coef(m_LBW_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue", xlab = "120 : 2016-12, 130: 2017-10, 140: 2018-08, 150: 2019-06, 160: 2020-04") #birth_Y_M_num
k.check(m_LBW_bam_covid_2)## k' edf k-index p-value
## s(GA_weeks) 9 7.292818 0.9667546 0.9375
## s(mat_age) 9 4.829053 0.9499728 0.5750
## s(birth_Y_M_num) 2 1.705205 0.9542470 0.7075 concurvity(m_LBW_bam_covid_2)## para s(GA_weeks) s(mat_age) s(birth_Y_M_num)
## worst 0.9828781 0.018060995 0.10773748 0.0011409949
## observed 0.9828781 0.009161863 0.04775956 0.0007068171
## estimate 0.9828781 0.009997521 0.06895348 0.0010619952 #to calculate the OR of LBW for non smooth terms
OR <- exp(coef(m_LBW_bam_covid_2))
beta <- coef(m_LBW_bam_covid_2)
Vb <- vcov(m_LBW_bam_covid_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 14), 2)## OR lci uci d
## (Intercept) 0.02 0.02 0.03 2.09
## mean_ssep2 0.99 0.99 0.99 0.00
## parity_cat(1,2] 0.50 0.48 0.53 0.38
## parity_cat(2,3] 0.42 0.39 0.46 0.48
## parity_cat(3,20] 0.45 0.39 0.51 0.44
## Urban3 1.02 0.98 1.07 0.01
## LanguageFrench 1.34 1.28 1.41 0.16
## LanguageItalian 1.18 1.05 1.32 0.09
## mother_nationality_cat2Africa 1.08 0.95 1.22 0.04
## mother_nationality_cat2Asia 0.88 0.78 0.98 0.07
## mother_nationality_cat2Europe 0.83 0.79 0.87 0.10
## mother_nationality_cat2Northern America 0.59 0.42 0.82 0.29
## mother_nationality_cat2Oceania 0.36 0.14 0.95 0.56
## mother_nationality_cat2Southern and Central America 0.77 0.65 0.91 0.15LBW vs “Normal” birthweight
2007-2020
Month was put as a linear variable.
m_LBW_normal_bam_2 <-bam(LBW_norm~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) +s(mean_ssep2)+ month + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7)
summary(m_LBW_normal_bam_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## LBW_norm ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) +
## month + parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.0369552 0.0196740 -205.192 < 2e-16 ***
## month 0.0007164 0.0017715 0.404 0.6859
## parity_cat(1,2] -0.6409693 0.0140031 -45.773 < 2e-16 ***
## parity_cat(2,3] -0.7685068 0.0227484 -33.783 < 2e-16 ***
## parity_cat(3,20] -0.7778748 0.0387627 -20.068 < 2e-16 ***
## Urban3 0.0181157 0.0126123 1.436 0.1509
## LanguageFrench 0.2804620 0.0140713 19.932 < 2e-16 ***
## LanguageItalian 0.1715774 0.0302897 5.665 1.47e-08 ***
## mother_nationality_cat2Africa 0.0115914 0.0377714 0.307 0.7589
## mother_nationality_cat2Asia -0.0143067 0.0317447 -0.451 0.6522
## mother_nationality_cat2Europe -0.1445716 0.0139627 -10.354 < 2e-16 ***
## mother_nationality_cat2Northern America -0.4899404 0.0912990 -5.366 8.04e-08 ***
## mother_nationality_cat2Oceania -0.4171044 0.2026645 -2.058 0.0396 *
## mother_nationality_cat2Southern and Central America -0.3583367 0.0458057 -7.823 5.16e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 8.134 8.636 86595.52 <2e-16 ***
## s(mat_age) 5.636 6.616 117.38 <2e-16 ***
## s(birth_Y_M_num) 3.911 4.836 54.43 <2e-16 ***
## s(mean_ssep2) 2.325 2.984 132.31 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.408 Deviance explained = 44.5%
## fREML = 1.4478e+06 Scale est. = 1 n = 1011626 plot(m_LBW_normal_bam_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_normal_bam_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks
plot(m_LBW_normal_bam_2, trans = plogis, select=2, ylim=c(.005, .03), shift = coef(m_LBW_normal_bam_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age
plot(m_LBW_normal_bam_2, trans = plogis, select=3, ylim=c(.005, .025), shift = coef(m_LBW_normal_bam_2)[1], shade = TRUE, shade.col = "lightblue", xlab = "1 : 2007-01, 25: 2009-01, 50: 2011-02, 75: 2013-03, 100: 2015-04, 125: 2017-05, 150: 2019-06") #birth_Y_M_num
plot(m_LBW_normal_bam_2, trans = plogis, select=4, ylim=c(.005, .030), shift = coef(m_LBW_normal_bam_2)[1], shade = TRUE, shade.col = "lightblue") #ssep
k.check(m_LBW_normal_bam_2)## k' edf k-index p-value
## s(GA_weeks) 9 8.133739 0.9273749 0.0225
## s(mat_age) 9 5.635822 0.9580069 0.4750
## s(birth_Y_M_num) 9 3.910756 0.9417167 0.1175
## s(mean_ssep2) 9 2.325108 0.9554838 0.3825 concurvity(m_LBW_normal_bam_2)## para s(GA_weeks) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst 0.8728924 0.018337743 0.12132443 0.03011043 0.17987362
## observed 0.8728924 0.007720419 0.04750016 0.01120391 0.12435460
## estimate 0.8728924 0.009454512 0.07649425 0.01162304 0.08214456 #to calculate the OR of LBW vs normal BW
OR <- exp(coef(m_LBW_normal_bam_2))
beta <- coef(m_LBW_normal_bam_2)
Vb <- vcov(m_LBW_normal_bam_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 14), 2)## OR lci uci d
## (Intercept) 0.02 0.02 0.02 2.23
## month 1.00 1.00 1.00 0.00
## parity_cat(1,2] 0.53 0.51 0.54 0.35
## parity_cat(2,3] 0.46 0.44 0.48 0.42
## parity_cat(3,20] 0.46 0.43 0.50 0.43
## Urban3 1.02 0.99 1.04 0.01
## LanguageFrench 1.32 1.29 1.36 0.15
## LanguageItalian 1.19 1.12 1.26 0.09
## mother_nationality_cat2Africa 1.01 0.94 1.09 0.01
## mother_nationality_cat2Asia 0.99 0.93 1.05 0.01
## mother_nationality_cat2Europe 0.87 0.84 0.89 0.08
## mother_nationality_cat2Northern America 0.61 0.51 0.73 0.27
## mother_nationality_cat2Oceania 0.66 0.44 0.98 0.23
## mother_nationality_cat2Southern and Central America 0.70 0.64 0.76 0.202008-2010
As always, month wasnt kept in the model (correlation between month and birth_Y_M_num variables). The variables birth_Y_M_num and ssep were set as linear terms.
m_LBW_normal_bam_2009_2 <-bam(LBW_norm~s(GA_weeks)+ s(mat_age) + birth_Y_M_num + mean_ssep2 + parity_cat +Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2008_10)
summary(m_LBW_normal_bam_2009_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## LBW_norm ~ s(GA_weeks) + s(mat_age) + birth_Y_M_num + mean_ssep2 +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.659351 0.109561 -33.400 < 2e-16 ***
## birth_Y_M_num 0.007519 0.001296 5.802 6.57e-09 ***
## mean_ssep2 -0.009721 0.001648 -5.898 3.68e-09 ***
## parity_cat(1,2] -0.613393 0.031052 -19.754 < 2e-16 ***
## parity_cat(2,3] -0.704519 0.050583 -13.928 < 2e-16 ***
## parity_cat(3,20] -0.745120 0.085821 -8.682 < 2e-16 ***
## Urban3 0.055127 0.027710 1.989 0.046656 *
## LanguageFrench 0.240118 0.031344 7.661 1.85e-14 ***
## LanguageItalian 0.294844 0.061125 4.824 1.41e-06 ***
## mother_nationality_cat2Africa -0.059372 0.095736 -0.620 0.535150
## mother_nationality_cat2Asia -0.024569 0.076027 -0.323 0.746571
## mother_nationality_cat2Europe -0.092698 0.031067 -2.984 0.002847 **
## mother_nationality_cat2Northern America -0.725093 0.208405 -3.479 0.000503 ***
## mother_nationality_cat2Oceania -0.138224 0.388947 -0.355 0.722305
## mother_nationality_cat2Southern and Central America -0.414559 0.099779 -4.155 3.26e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 7.077 7.810 17825.65 < 2e-16 ***
## s(mat_age) 3.840 4.752 21.52 0.000506 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.4 Deviance explained = 43.3%
## fREML = 2.8729e+05 Scale est. = 1 n = 201968 plot(m_LBW_normal_bam_2009_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_normal_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks
plot(m_LBW_normal_bam_2009_2, trans = plogis, select=2, ylim=c(.005, .05), shift = coef(m_LBW_normal_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age
k.check(m_LBW_normal_bam_2009_2)## k' edf k-index p-value
## s(GA_weeks) 9 7.077019 0.9182415 0.010
## s(mat_age) 9 3.839971 0.9712118 0.965 concurvity(m_LBW_normal_bam_2009_2)## para s(GA_weeks) s(mat_age)
## worst 0.9847922 0.016590747 0.13530354
## observed 0.9847922 0.006969109 0.04023792
## estimate 0.9847922 0.008342767 0.08762967 #to calculate the OR of LBW vs normal BW
OR <- exp(coef(m_LBW_normal_bam_2009_2))
beta <- coef(m_LBW_normal_bam_2009_2)
Vb <- vcov(m_LBW_normal_bam_2009_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 15), 2)## OR lci uci d
## (Intercept) 0.03 0.02 0.03 2.02
## birth_Y_M_num 1.01 1.00 1.01 0.00
## mean_ssep2 0.99 0.99 0.99 0.01
## parity_cat(1,2] 0.54 0.51 0.58 0.34
## parity_cat(2,3] 0.49 0.45 0.55 0.39
## parity_cat(3,20] 0.47 0.40 0.56 0.41
## Urban3 1.06 1.00 1.12 0.03
## LanguageFrench 1.27 1.20 1.35 0.13
## LanguageItalian 1.34 1.19 1.51 0.16
## mother_nationality_cat2Africa 0.94 0.78 1.14 0.03
## mother_nationality_cat2Asia 0.98 0.84 1.13 0.01
## mother_nationality_cat2Europe 0.91 0.86 0.97 0.05
## mother_nationality_cat2Northern America 0.48 0.32 0.73 0.40
## mother_nationality_cat2Oceania 0.87 0.41 1.87 0.08
## mother_nationality_cat2Southern and Central America 0.66 0.54 0.80 0.232017-2020
Ssep variable was set as linear variables (when included as smooth variable, edf= 1.62)
m_LBW_normal_bam_covid_2 <-bam(LBW_norm~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) + mean_ssep2 + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2017_20)
summary(m_LBW_normal_bam_covid_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## LBW_norm ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + mean_ssep2 +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.618770 0.085835 -42.160 < 2e-16 ***
## mean_ssep2 -0.008824 0.001379 -6.400 1.55e-10 ***
## parity_cat(1,2] -0.675156 0.025787 -26.182 < 2e-16 ***
## parity_cat(2,3] -0.840120 0.041915 -20.044 < 2e-16 ***
## parity_cat(3,20] -0.780183 0.068901 -11.323 < 2e-16 ***
## Urban3 0.022726 0.023047 0.986 0.32410
## LanguageFrench 0.290190 0.025667 11.306 < 2e-16 ***
## LanguageItalian 0.153307 0.059736 2.566 0.01028 *
## mother_nationality_cat2Africa 0.081889 0.063870 1.282 0.19980
## mother_nationality_cat2Asia -0.131873 0.057235 -2.304 0.02122 *
## mother_nationality_cat2Europe -0.182846 0.025529 -7.162 7.94e-13 ***
## mother_nationality_cat2Northern America -0.524457 0.173262 -3.027 0.00247 **
## mother_nationality_cat2Oceania -1.026005 0.489925 -2.094 0.03624 *
## mother_nationality_cat2Southern and Central America -0.256735 0.085221 -3.013 0.00259 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 7.269 7.907 25423.951 <2e-16 ***
## s(mat_age) 4.835 5.826 55.708 <2e-16 ***
## s(birth_Y_M_num) 1.856 2.316 3.334 0.262
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.409 Deviance explained = 45.4%
## fREML = 4.3553e+05 Scale est. = 1 n = 306277 plot(m_LBW_normal_bam_covid_2, trans = plogis, select=1, ylim=c(0, 1), shift = coef(m_LBW_normal_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks
plot(m_LBW_normal_bam_covid_2, trans = plogis, select=2, ylim=c(.005, .07), shift = coef(m_LBW_normal_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue") #mat_age
plot(m_LBW_normal_bam_covid_2, trans = plogis, select=3, ylim=c(.020, .028), shift = coef(m_LBW_normal_bam_covid_2)[1], shade = TRUE, shade.col = "lightblue", xlab = "120 : 2016-12, 130: 2017-10, 140: 2018-08, 150: 2019-06, 160: 2020-04") #birth_Y_M_num
k.check(m_LBW_normal_bam_covid_2)## k' edf k-index p-value
## s(GA_weeks) 9 7.269236 0.9003570 0.0000
## s(mat_age) 9 4.835330 0.9640484 0.9425
## s(birth_Y_M_num) 9 1.855514 0.9436354 0.4875 concurvity(m_LBW_normal_bam_covid_2)## para s(GA_weeks) s(mat_age) s(birth_Y_M_num)
## worst 0.9827989 0.020008547 0.10645213 0.001371803
## observed 0.9827989 0.009105054 0.04684491 0.001008478
## estimate 0.9827989 0.010278237 0.06782008 0.001042187#to calculate the OR of LBW vs normal BW
OR <- exp(coef(m_LBW_normal_bam_covid_2))
beta <- coef(m_LBW_normal_bam_covid_2)
Vb <- vcov(m_LBW_normal_bam_covid_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 14), 2)## OR lci uci d
## (Intercept) 0.03 0.02 0.03 2.00
## mean_ssep2 0.99 0.99 0.99 0.00
## parity_cat(1,2] 0.51 0.48 0.54 0.37
## parity_cat(2,3] 0.43 0.40 0.47 0.46
## parity_cat(3,20] 0.46 0.40 0.52 0.43
## Urban3 1.02 0.98 1.07 0.01
## LanguageFrench 1.34 1.27 1.41 0.16
## LanguageItalian 1.17 1.04 1.31 0.08
## mother_nationality_cat2Africa 1.09 0.96 1.23 0.05
## mother_nationality_cat2Asia 0.88 0.78 0.98 0.07
## mother_nationality_cat2Europe 0.83 0.79 0.88 0.10
## mother_nationality_cat2Northern America 0.59 0.42 0.83 0.29
## mother_nationality_cat2Oceania 0.36 0.14 0.94 0.57
## mother_nationality_cat2Southern and Central America 0.77 0.65 0.91 0.14Macrosomia vs “Normal” birthweight
2007-2020
m_macro_normal_bam <-bam(macro_norm~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) + s(month) +s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7)
summary(m_macro_normal_bam)##
## Family: binomial
## Link function: logit
##
## Formula:
## macro_norm ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + s(month) +
## s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.988756 0.008479 -352.494 < 2e-16 ***
## parity_cat(1,2] 0.582219 0.008303 70.125 < 2e-16 ***
## parity_cat(2,3] 0.841384 0.011433 73.594 < 2e-16 ***
## parity_cat(3,20] 0.987128 0.018273 54.020 < 2e-16 ***
## Urban3 0.013847 0.007668 1.806 0.07095 .
## LanguageFrench -0.268400 0.009394 -28.571 < 2e-16 ***
## LanguageItalian -0.372239 0.023024 -16.168 < 2e-16 ***
## mother_nationality_cat2Africa 0.061706 0.021283 2.899 0.00374 **
## mother_nationality_cat2Asia -0.109246 0.022227 -4.915 8.88e-07 ***
## mother_nationality_cat2Europe 0.171988 0.008118 21.186 < 2e-16 ***
## mother_nationality_cat2Northern America 0.337389 0.043411 7.772 7.73e-15 ***
## mother_nationality_cat2Oceania 0.181090 0.104897 1.726 0.08428 .
## mother_nationality_cat2Southern and Central America 0.148646 0.028587 5.200 1.99e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 7.728 8.431 30433.750 <2e-16 ***
## s(mat_age) 4.525 5.449 12.595 0.0369 *
## s(birth_Y_M_num) 4.301 5.295 77.959 <2e-16 ***
## s(month) 2.842 3.532 8.933 0.0468 *
## s(mean_ssep2) 6.664 7.813 61.399 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.0445 Deviance explained = 7.43%
## fREML = 1.4888e+06 Scale est. = 1 n = 1050380plot(m_macro_normal_bam, trans = plogis, select=1, ylim=c(0, 0.6), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks
plot(m_macro_normal_bam, trans = plogis, select=2, ylim=c(0.03, 0.06), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue") #mat_age
plot(m_macro_normal_bam, trans = plogis, select=3, ylim=c(0.045, 0.055), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue", xlab = "1 : 2007-01, 25: 2009-01, 50: 2011-02, 75: 2013-03, 100: 2015-04, 125: 2017-05, 150: 2019-06") #birth_Y_M_num
plot(m_macro_normal_bam, trans = plogis, select=4, ylim=c(0.045, 0.052), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue") #month
plot(m_macro_normal_bam, trans = plogis, select=5, ylim=c(0.02, 0.06), shift = coef(m_macro_normal_bam)[1], shade = TRUE, shade.col = "lightblue") #ssep
k.check(m_macro_normal_bam)## k' edf k-index p-value
## s(GA_weeks) 9 7.728394 0.9504278 0.8575
## s(mat_age) 9 4.525169 0.9498593 0.8175
## s(birth_Y_M_num) 9 4.301243 0.9400573 0.5875
## s(month) 9 2.842075 0.9521674 0.8775
## s(mean_ssep2) 9 6.664473 0.9509826 0.8350 concurvity(m_macro_normal_bam)## para s(GA_weeks) s(mat_age) s(birth_Y_M_num) s(month) s(mean_ssep2)
## worst 0.7583176 0.01888534 0.125009429 0.030368293 0.0282971647 0.17894841
## observed 0.7583176 0.01351314 0.007078764 0.009518221 0.0008605803 0.02650765
## estimate 0.7583176 0.01216139 0.080052791 0.011622466 0.0211814067 0.08691995#to calculate the OR of Macrosomia vs normal BW
OR <- exp(coef(m_macro_normal_bam))
beta <- coef(m_macro_normal_bam)
Vb <- vcov(m_macro_normal_bam, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 13), 2)## OR lci uci d
## (Intercept) 0.05 0.05 0.05 1.65
## parity_cat(1,2] 1.79 1.76 1.82 0.32
## parity_cat(2,3] 2.32 2.27 2.37 0.46
## parity_cat(3,20] 2.68 2.59 2.78 0.55
## Urban3 1.01 1.00 1.03 0.01
## LanguageFrench 0.76 0.75 0.78 0.15
## LanguageItalian 0.69 0.66 0.72 0.21
## mother_nationality_cat2Africa 1.06 1.02 1.11 0.03
## mother_nationality_cat2Asia 0.90 0.86 0.94 0.06
## mother_nationality_cat2Europe 1.19 1.17 1.21 0.10
## mother_nationality_cat2Northern America 1.40 1.29 1.53 0.19
## mother_nationality_cat2Oceania 1.20 0.98 1.47 0.10
## mother_nationality_cat2Southern and Central America 1.16 1.10 1.23 0.082008-2010
We set mat age and birth_Y_M_num variables as linear (edf= 1.005 and 1.378, respectively)
m_macro_normal_bam_2009_2 <-bam(macro_norm~s(GA_weeks)+ s(mean_ssep2) + mat_age + birth_Y_M_num +parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2008_10)
summary(m_macro_normal_bam_2009_2)##
## Family: binomial
## Link function: logit
##
## Formula:
## macro_norm ~ s(GA_weeks) + s(mean_ssep2) + mat_age + birth_Y_M_num +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.0624066 0.0594550 -51.508 < 2e-16 ***
## mat_age 0.0033146 0.0017206 1.926 0.05406 .
## birth_Y_M_num -0.0008767 0.0007783 -1.126 0.25997
## parity_cat(1,2] 0.5519526 0.0184340 29.942 < 2e-16 ***
## parity_cat(2,3] 0.8393671 0.0252785 33.205 < 2e-16 ***
## parity_cat(3,20] 0.9752223 0.0400244 24.366 < 2e-16 ***
## Urban3 0.0148538 0.0168953 0.879 0.37931
## LanguageFrench -0.2539002 0.0210518 -12.061 < 2e-16 ***
## LanguageItalian -0.3262404 0.0478543 -6.817 9.27e-12 ***
## mother_nationality_cat2Africa 0.3452120 0.0536319 6.437 1.22e-10 ***
## mother_nationality_cat2Asia -0.1520186 0.0536794 -2.832 0.00463 **
## mother_nationality_cat2Europe 0.1923115 0.0182076 10.562 < 2e-16 ***
## mother_nationality_cat2Northern America 0.2987413 0.0954085 3.131 0.00174 **
## mother_nationality_cat2Oceania 0.1256268 0.2318870 0.542 0.58798
## mother_nationality_cat2Southern and Central America 0.2611198 0.0577739 4.520 6.19e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 6.602 7.499 6418.525 <2e-16 ***
## s(mean_ssep2) 2.339 3.010 6.209 0.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.047 Deviance explained = 7.81%
## fREML = 2.9711e+05 Scale est. = 1 n = 209945 plot(m_macro_normal_bam_2009_2, trans = plogis, select=1, ylim=c(0, 0.8), shift = coef(m_macro_normal_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks
plot(m_macro_normal_bam_2009_2, trans = plogis, select=2, ylim=c(0.03, 0.06), shift = coef(m_macro_normal_bam_2009_2)[1], shade = TRUE, shade.col = "lightblue") #ssep
k.check(m_macro_normal_bam_2009_2)## k' edf k-index p-value
## s(GA_weeks) 9 6.60227 0.9345813 0.4100
## s(mean_ssep2) 9 2.33933 0.9479809 0.7875 concurvity(m_macro_normal_bam_2009_2)## para s(GA_weeks) s(mean_ssep2)
## worst 0.9806332 0.01626321 0.17850048
## observed 0.9806332 0.01135410 0.16053032
## estimate 0.9806332 0.01127897 0.08817556#to calculate the OR of Macrosomia vs normal BW
OR <- exp(coef(m_macro_normal_bam_2009_2))
beta <- coef(m_macro_normal_bam_2009_2)
Vb <- vcov(m_macro_normal_bam_2009_2, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 15), 2)## OR lci uci d
## (Intercept) 0.05 0.04 0.05 1.69
## mat_age 1.00 1.00 1.01 0.00
## birth_Y_M_num 1.00 1.00 1.00 0.00
## parity_cat(1,2] 1.74 1.68 1.80 0.30
## parity_cat(2,3] 2.31 2.20 2.43 0.46
## parity_cat(3,20] 2.65 2.45 2.87 0.54
## Urban3 1.01 0.98 1.05 0.01
## LanguageFrench 0.78 0.74 0.81 0.14
## LanguageItalian 0.72 0.66 0.79 0.18
## mother_nationality_cat2Africa 1.41 1.27 1.57 0.19
## mother_nationality_cat2Asia 0.86 0.77 0.95 0.08
## mother_nationality_cat2Europe 1.21 1.17 1.26 0.11
## mother_nationality_cat2Northern America 1.35 1.12 1.63 0.17
## mother_nationality_cat2Oceania 1.13 0.72 1.79 0.07
## mother_nationality_cat2Southern and Central America 1.30 1.16 1.45 0.142017-2020
m_macro_normal_bam_covid <-bam(macro_norm~s(GA_weeks)+ s(mat_age) + s(birth_Y_M_num) +s(mean_ssep2) + parity_cat + Urban3 + Language + mother_nationality_cat2,family="binomial", data=bevn_eco_in7_2017_20)
summary(m_macro_normal_bam_covid)##
## Family: binomial
## Link function: logit
##
## Formula:
## macro_norm ~ s(GA_weeks) + s(mat_age) + s(birth_Y_M_num) + s(mean_ssep2) +
## parity_cat + Urban3 + Language + mother_nationality_cat2
##
## Parametric coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.966355 0.015301 -193.867 < 2e-16 ***
## parity_cat(1,2] 0.614600 0.015109 40.677 < 2e-16 ***
## parity_cat(2,3] 0.875458 0.020714 42.264 < 2e-16 ***
## parity_cat(3,20] 1.019343 0.033436 30.486 < 2e-16 ***
## Urban3 0.005643 0.013881 0.407 0.684368
## LanguageFrench -0.289867 0.016778 -17.277 < 2e-16 ***
## LanguageItalian -0.421656 0.045138 -9.342 < 2e-16 ***
## mother_nationality_cat2Africa -0.051770 0.035678 -1.451 0.146767
## mother_nationality_cat2Asia -0.105231 0.038140 -2.759 0.005796 **
## mother_nationality_cat2Europe 0.130878 0.014712 8.896 < 2e-16 ***
## mother_nationality_cat2Northern America 0.299535 0.083080 3.605 0.000312 ***
## mother_nationality_cat2Oceania -0.058639 0.242115 -0.242 0.808630
## mother_nationality_cat2Southern and Central America 0.116127 0.056445 2.057 0.039651 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df Chi.sq p-value
## s(GA_weeks) 6.203 7.075 8864.53 < 2e-16 ***
## s(mat_age) 3.889 4.782 18.33 0.00175 **
## s(birth_Y_M_num) 2.444 3.044 12.81 0.00524 **
## s(mean_ssep2) 4.463 5.563 41.08 1.85e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.043 Deviance explained = 7.15%
## fREML = 4.5124e+05 Scale est. = 1 n = 318299plot(m_macro_normal_bam_covid, trans = plogis, select=1, ylim=c(0, 0.7), shift = coef(m_macro_normal_bam_covid)[1], shade = TRUE, shade.col = "lightblue") #GA_weeks
plot(m_macro_normal_bam_covid, trans = plogis, select=2, ylim=c(0.02, 0.08), shift = coef(m_macro_normal_bam_covid)[1], shade = TRUE, shade.col = "lightblue") #mat_age
plot(m_macro_normal_bam_covid, trans = plogis, select=3, ylim=c(0.04, 0.06), shift = coef(m_macro_normal_bam_covid)[1], shade = TRUE, shade.col = "lightblue", xlab = "120 : 2016-12, 130: 2017-10, 140: 2018-08, 150: 2019-06, 160: 2020-04") #birth_Y_M_num
plot(m_macro_normal_bam_covid, trans = plogis, select=4, ylim=c(0.02, 0.06), shift = coef(m_macro_normal_bam_covid)[1], shade = TRUE, shade.col = "lightblue") #ssep
k.check(m_macro_normal_bam_covid)## k' edf k-index p-value
## s(GA_weeks) 9 6.203012 0.9327167 0.1775
## s(mat_age) 9 3.888923 0.9544342 0.7525
## s(birth_Y_M_num) 9 2.444071 0.9590803 0.8325
## s(mean_ssep2) 9 4.463401 0.9251614 0.0475 concurvity(m_macro_normal_bam_covid)## para s(GA_weeks) s(mat_age) s(birth_Y_M_num) s(mean_ssep2)
## worst 0.7608652 0.02085818 0.11064626 0.0013996579 0.17772643
## observed 0.7608652 0.01556035 0.07319894 0.0008908752 0.01327785
## estimate 0.7608652 0.01460350 0.07032166 0.0009817982 0.07889043#to calculate the OR of Macrosomia vs normal BW
OR <- exp(coef(m_macro_normal_bam_covid))
beta <- coef(m_macro_normal_bam_covid)
Vb <- vcov(m_macro_normal_bam_covid, unconditional = TRUE)
se <- sqrt(diag(Vb))
lci <- exp(beta-1.96*se)
uci <- exp(beta+1.96*se)
d <- abs(beta)/1.81
my.ci <- data.frame(cbind(OR, lci, uci, d))
round(head(my.ci, 13), 2)## OR lci uci d
## (Intercept) 0.05 0.05 0.05 1.64
## parity_cat(1,2] 1.85 1.79 1.90 0.34
## parity_cat(2,3] 2.40 2.30 2.50 0.48
## parity_cat(3,20] 2.77 2.60 2.96 0.56
## Urban3 1.01 0.98 1.03 0.00
## LanguageFrench 0.75 0.72 0.77 0.16
## LanguageItalian 0.66 0.60 0.72 0.23
## mother_nationality_cat2Africa 0.95 0.89 1.02 0.03
## mother_nationality_cat2Asia 0.90 0.84 0.97 0.06
## mother_nationality_cat2Europe 1.14 1.11 1.17 0.07
## mother_nationality_cat2Northern America 1.35 1.15 1.59 0.17
## mother_nationality_cat2Oceania 0.94 0.59 1.52 0.03
## mother_nationality_cat2Southern and Central America 1.12 1.01 1.25 0.06